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SUMMARY 

An  exponential  finite  difference  scheme  first  presented  by  Bhattacharya 
for  one-dimensional  unsteady  heat  conduction  problems  in  Cartesian  coordinates 
has  been  extended.  The  finite  difference  algorithm  developed  was  used  to  solve 
the  unsteady  diffusion  equation  In  one-dimensional  cylindrical  coordinates  and 
was  applied  to  two-  and  three-dimensional  conduction  problems  in  Cartesian 
coordinates.  Heat  conduction  involving  variable  thermal  conductivity  was  also 
investigated.  The  method  was  used  to  solve  nonlinear  partial  differential 
equations  in  one-  (Burger's  equation)  and  two-  (boundary  layer  equations) 
dimensional  Cartesian  coordinates.  Predicted  results  are  compared  to  exact 
solutions  where  available  or  to  results  obtained  by  other  numerical  methods. 


INTRODUCTION 

The  objective  of  this  work  is  to  extend,  expand,  and  compare  an  explicit 
exponential  finite  difference  technique  first  proposed  by  Bhattacharya 
(ref.  1).  To  date  the  method  has  only  been  used  for  one-dimensional  unsteady 
heat  transfer  in  Cartesian  coordinates.  The  method  is  a  finite  difference  rel¬ 
ative  of  the  separation  of  variables  technique.  The  finite  difference  equa¬ 
tion  that  results  uses  time  step  division  to  increase  accuracy  and  to  maintain 
stabi 1 i ty. 

Following  his  initial  paper,  Bhattacharya  (ref.  2)  and  Bhattacharya  and 
Davies  (ref.  3)  have  developed  refined  forms  of  the  exponential  finite  differ¬ 
ence  equation.  Also,  an  approximate  substitution  for  a  given  range  of  expo¬ 
nential  term  was  Investigated  to  reduce  the  computation  time  while  retaining 
good  accuracy.  In  references  1  to  3,  the  results  for  unsteady  one-dimensional 
heat  transfer  found  by  implicit  and  explicit  numerical  techniques  were  com¬ 
pared  to  exact  analysis.  The  overall  results  indicated  that  the  exponential 
finite  difference  techniques  were  more  accurate  than  the  other  available  numer¬ 
ical  techniques.  The  one  drawback  with  the  exponential  finite  difference 
method  was  that  computer  time  increased  for  the  one-dimensional  case  that  was 
investigated. 

The  intent  of  the  present  work  is  to  demonstrate  how  the  exponential 
finite  difference  method  originally  developed  in  reference  1  can  be  used  to 
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solve  a  wide  variety  of  problems.  Linear  and  nonlinear  partial  differential 
equations  found  In  engineering  and  physics  will  be  solved.  All  results 
obtained  by  this  finite  difference  technique  will  be  compared  to  exact  solu¬ 
tions  or  to  values  found  by  use  of  other  numerical  techniques. 
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NOMENCLATURE 


Biot  modulus,  hR/k 


material  specific  heat,  J / ( kg > < <> ;  Btu/( 1 bm)(°F) 

convection  heat  transfer  coefficient,  W/(m2)(°C);  Btu/(ft2)(hr)(°F) 

nodal  location  in  x,  y,  and  z  spatial  coordinate  directions, 
respectively 

Bessel  functions  of  zero  and  first  order,  respectively 
thermal  conductivity,  W/(m)(°C);  Btu/(hr)(°F)(ft) 

thermal  conductivity  at  1th  position,  nth  time  step,  W/(m)(°C); 
Btu/(hr)(°F)(ft) 

distance  between  plates,  m;  ft 

dimensionless  drive  number 

number  of  subintervals 

number  of  nodes  In  a  spatial  direction 

time  step  position  designation 


R,Rl,R2  radial  length,  m;  ft 


spatial  coordinate  (cylindrical  coordinates),  m;  ft 
temperature,  °C;  °F 
time,  s 

time  between  time  steps  n  and  n  +  1 

flow  velocity,  m/s;  ft/s 

method  of  Douglas  Intermediate  values 

substitution  variable  for  Burger's  equation 

spatial  coordinates  (Cartesian  coordinates),  m;  ft 

distance  between  nodal  positions  in  the  x,  y,  and  z  spatial 
directions,  respectively 
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thermal  dlffuslvlty,  m2/s;  ft2/s 
rate  of  thermal  conductivity  variation 
At/pCp(Ax)2,  [W/(m)(°C)]"1 ;  CBtu/(hr) (°f )(ft2) ]_1 
constant 

constant  used  in  exponential  finite  difference  method  with  tempera¬ 
ture-varying  thermal  conductivity 

finite  difference  operator 

mt*1  eigenvalue  of  Bessel  function 

amplification  factor 

kinematic  viscosity,  m2/s;  ft2 / s 

material  density,  kg/m^;  lbm/ft^ 

Kirchoff  transformation  variable 

(a  At)/(Ax)2  dimensionless  time 

separation  variables 


ANALYSIS 

The  exponential  finite  difference  algorithm  derived  by  Bhattacharya 
(ref.  1)  will  be  developed  in  this  section.  To  illustrate  the  procedure 
unsteady  two-dimensional  heat  conduction  in  Cartesian  coordinates  will  be 
initially  considered  (ref.  4).  The  appropriate  partial  differential  equation 
is  (ref.  5): 
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where  a  is  the  thermal  diffusivity,  k/pc.  If  equation  (1)  is  divided  by  T 
and  the  resulting  expression  evaluated  at  the  time  step  n  and  the  grid  point 
( i , j ) ,  we  may  write 


9  In  T 

at 


a  /Vt  3 h 

T  lax2  +  ay2 


It  may  be  assumed  here  that  T  can  be  written  in  a  product  form  as 

T  =  $(x,y)  g/< t ) 
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The  spatial  and  temporal  variables  can  be  separated  and  equation  (2)  can  then 
be  set  equal  to  a  constant,  say,  -k.  Thus,  the  left  side  of  equation  (2)  Is 
written  as 

3  In  T  n 

Replacing  the  derivative  with  a  one-sided  difference  results  in  the  following: 
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The  separation  constant  k  Is  evaluated  from  the  right  side  of  equation  (2) 
by  using  second  central  space  differences  which  result  In 
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If  the  grid  spacing  Is  constant  <Ax  =  Ay),  then  equation  (3)  may  be  written  as 
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Q  =  grid  Fourier  number  =  (a  At)/Ax  and  M,  . 
r  which  is  written  as 


=  dimensionless  drive 
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Because  of  the  exponential  form  of  equation  (5),  the  time  step  may  be 
divided  into  a  number  of  subintervals.  Subdivisions  or  reduction  of  the  time 
step  is  typically  done  to  increase  the  accuracy  of  explicit  numerical  methods. 

For  example,  if  the  time  step  were  divided  into  two  intervals,  then  T?+l 
would  be  found  In  the  following  way:  ’ ^ 
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where  the  dimensionless  drive  numbers  are  evaluated  at  the  sub-time  intervals 
and  then  summed  for  calculation  of  "T"  at  the  n  +  1  time  step.  Or  in  a  more 
general  form,  for  "m"  subintervals 
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(7) 


Equation  (7)  is  the  general  difference  equation  for  the  temperature  at  the  1,j 
node,  at  the  n+1  time  step,  for  m  time-step  subintervals.  This  equation 
is  valid  for  all  interior  nodes  for  two-dimensional  rectangular  domain.  Nodes 
on  the  boundaries  are  treated  differently  and  depend  on  the  application. 


In  reference  1,  it  was  shown  that  for  heat  transfer  applications  the  time 
step  can  be  subdivided  to  a  maximum  number  of  time  subintervals  as  follows: 


S(N/2)  -  1  ♦  heat  transfer  coefficient  »  infinite 

(8) 

(N/2)  +  1  ♦  heat  transfer  coefficient  finite 
where  N  equals  the  number  of  nodes  in  one  of  the  coordinate  directions. 


STABILITY  OF  THE  EXPONENTIAL  FINITE  DIFFERENCE  METHOD 

With  few  exceptions,  explicit  finite  difference  procedures  for  solving 
partial  differential  equations  are  inherently  unstable  unless  certain  numerical 
conditions  are  satisfied.  These  conditions  take  the  form  of  a  grid  size  and/or 
time  step  requirement  written  in  terms  of  parameters  of  the  given  problem.  If 
these  stability  conditions  are  not  met,  the  solution  can  diverge.  These  con¬ 
straints  on  grid  size  or  length  of  time  step  can  make  the  methods  impractical 
for  certain  applications.  These  conditions,  however,  must  be  known  prior  to 
use  of  any  explicit  differencing  procedure. 

There  are  a  variety  of  methods  that  have  been  used  to  establish  tl.e  sta¬ 
bility  constraints  of  a  finite  difference  procedure.  These  methods  seek  to 
find  an  expression  for  the  amplification  factor  which  is  defined  as  the  ratio 
of  the  current  solution  result  to  that  in  the  previous  step.  If  the  absolute 
value  of  the  ratio  is  less  than  one,  then  the  method  is  regarded  as  being  sta¬ 
ble.  Determination  of  the  amplification  factor  for  the  exponential  finite  dif¬ 
ference  method  is  particularly  convenient,  as  has  been  shown  in  reference  1. 

For  the  two-dimensional  Cartesian  coordinate  case,  the  amplification  factor 
E  can  be  defined  as  the  following  (no  time  subinterval  division): 
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Numerical  stability  constraints  require  that 


11m  |d  <  1  (10) 

At-»0 

Ax-»0 

To  satisfy  this  requirement,  the  exponent  of  equation  (9)  must  obviously  be 
less  than  or  equal  to  zero.  Since  the  components  that  make  up  Q  in  that 
exponent  are  all  positive,  this  implies  that  the  dimensionless  drive  number 
will  determine  the  numerical  stability.  For  the  two-dimensional  Cartesian 
coordinate  case  the  dimensionless  drive  number  must  satisfy 
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Equation  (12)  needs  to  be  satisfied  otherwise  an  unstable  condition  can 
exist.  As  pointed  out  by  Bhattacharya  (ref.  1),  the  dimensionless  drive  number 
primarily  determines  the  stability  of  the  solution.  However  a  large  dimension¬ 
less  time  step  could  also  cause  the  solution  to  become  unstable.  Since  time 
subinterval  division  is  used,  the  total  dimensionless  time  step  Q  could 
become  quite  large.  In  reference  1,  it  was  recommended  for  one-dimensional 
heat  conduction  problems  that  the  dimensionless  time  step  satisfy  the  follow¬ 
ing  condition: 

~~t  <  0.5  (13) 


where  m  is  the  number  of  time-step  subintervals  involved  in  the  calculations. 
This  same  reasoning  can  be  extended  to  heat  conduction  problems  in  two  and 
three  dimensions  with  equal  grid  spacing.  The  expression  in  equation  (13)  has 
been  shown  in  reference  4  to  be  equal  to  1/4  and  1/6  for  two  and  three  dimen¬ 
sions,  respectively.  This  restriction  as  shown  in  equation  (13)  is  of  the  same 
magnitude  as  is  typically  used  for  an  explicit  finite  difference  technique  for 
the  grid  Fourier  number. 


APPLICATIONS 

The  exponential  finite  difference  technique  will  now  be  applied  to  a  num¬ 
ber  of  engineering  problems.  Unsteady  heat  transfer  problems  will  be  solved 
in  one-dimensional  radial  coordinates,  in  one-dimensional  Cartesian  coordi¬ 
nates  with  temperature-varying  thermal  conductivity,  and  in  three-dimensional 
Cartesian  coordinates.  Nonlinear  equations  will  also  be  numerically  solved 
using  this  method.  In  particular,  Burger's  equation  and  the  laminar  boundary 
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layer  on  a  flat  plate  will  be  Investigated.  All  applications  will  be  compared 
either  to  exact  results  or  to  results  obtained  via  other  numerical  techniques. 
This  comparison  will  provide  an  assessment  of  the  accuracy  of  the  exponential 
finite  difference  method. 


One-Dimensional  Heat  Conduction  In  Cylindrical  Coordinates 

One-dimensional  heat  conduction  In  cylindrical  coordinates  will  be  inves¬ 
tigated  for  Infinite  and  finite  heat  transfer  coefficient.  The  exact  results 
for  both  cases  can  be  found  In  reference  5. 

For  infinite  heat  transfer  coefficient  on  the  boundary  surface  the  exact 
result  Is  given  In  reference  5  as 
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where  X-R  Is  the  m^h  zero  of 
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Jn(X-R)  =  0 

o  m 


(15) 


The  results  of  both  the  exact  analysis  and  the  exponential  finite  differ¬ 
ence  method  are  shown  In  table  I.  As  can  be  seen  from  the  tabulated  results, 
exponential  finite  difference  results  approach  the  exact  solution  as  the  number 
of  nodes  Is  Increased  or  as  the  dimensionless  time  step  is  decreased. 

Hhen  the  heat  transfer  coefficient  has  a  finite  value  at  the  surface,  the 
exact  solution  from  reference  5  is 
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where  B  =  hR/k  (Biot  modulus)  and  X-  (characteristics  values)  are  given  by 

m  3  J 

the  following  equation  (for  cooling): 


(X-R) J, (X-R)  -  BJn(X-R)  =  0 
m  I  m  0  m 


(17) 


The  results  are  shown  in  table  II  for  various  values  of  the  Biot  modulus. 
As  would  be  expected,  the  solution  approaches  the  exact  solution  as  the  number 
of  nodes  Increases.  As  the  elapsed  time  of  the  solution  proceeded,  tempera¬ 
tures  predicted  by  the  exponential  finite  difference  method  approached  the 
exact  result.  Also  the  results  indicated  that  reducing  the  size  of  the  time 
subinterval  Increased  the  accuracy  of  the  method. 


One  last  comparison  will  be  made  while  Investigating  the  exponential 
finite  difference  technique  in  one-dimensional  cylindrical  coordinates.  The 
geometry  for  a  cylindrical  annulus  Is  shown  in  figure  1  and  is  applied  to  a 
problem  with  the  following  initial  and  boundary  conditions: 

T(r,0)  =  0 
T<R2,t>  =  1 .0 

If  (IV‘>  ■  0 

In  reference  6  this  problem  was  solved  numerically  using  a  character i sti c- 
value  solution.  A  comparison  of  results  is  shown  in  table  III  for  the  exponen¬ 
tial  method  using  the  same  grid  spacing  as  in  reference  6  and  for  the  case 
where  grid  spacing  is  halved.  The  results  are  seen  to  compare  quite  well  with 
the  finer  mesh  being  slightly  closer  to  the  value  from  reference  6  especially 
during  the  first  few  time  steps  of  the  solution. 


One-Dimensional  Unsteady  State  Conduction  With 
Temperature-Varying  Thermal  Conductivity 

The  effect  of  temperature-varying  thermal  conductivity  will  now  be  inves¬ 
tigated  using  three  different  numerical  schemes:  a  pure  explicit,  the  expo¬ 
nential  method,  and  an  implicit  technique.  The  problem  to  be  solved  is 
illustrated  in  figure  2(a).  The  thermal  conductivity  is  assumed  to  be  a  linear 
function  of  temperature  and  is  shown  in  figure  2(b). 

The  exponential  finite  difference  method  will  be  applied  first  to  the 
given  problem.  The  following  governing  partial  differential  equation  is  taken 
from  reference  7: 

”cp  If  -  If  (k  U)  (l9> 

Equation  (19)  can  be  changed  to  a  simpler  form  by  using  an  alternate 
dependent  variable  0  (the  Kirchoff  transformation)  given  by 
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Since  it  has  been  assumed  that  the  thermal  conductivity  is  a  linear  function 
of  temperature, 

k(T)  =  kR(l  +  (5T)  (2 

Now  substituting  equation  (23)  into  equation  (20)  results  in  the  following: 


0  -  jj-  (kR  +  |3kRT)dT 
R  Jj 
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Direct  integration  yields: 
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Equation  (24)  provides  the  relationship  between  the  variable  T  and  the  vari¬ 
able  0. 

Returning  to  equation  (22)  and  rearranging  results  in: 

ae_  _k_  afe_  (0, 
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Equation  (25)  is  in  a  form  for  which  the  exponential  finite  difference 
method  can  be  applied.  The  resulting  equation  in  the  Kirchoff  variable  can  be 
shown  (ref.  4)  to  be  given  by 
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Evaluating  equation  (24)  at  node  i  and  time  step  n  results  in 
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Substitution  of  equation  (27)  into  equation  (26)  at  the  appropriate  time  steps 
and  nodal  locations  yields 
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If  Tr  =  =  O.O,  equation  (28)  becomes 
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The  equation  for  T"  Is  a  quadratic  with  the  right  side  of  the  equation 
that  is  known  at  time  step  n.  Hence,  define  a  variable,  such  that 
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Equation  (29)  then  becomes 
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Solving  this  and  using  the  positive  root  results  in 
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where  (3  >  0. 


Equation  (32)  and  equation  (30)  are  solved  using  the  exponential  finite 
difference  solution  sequence.  In  this  case  the  conductivity  as  well  as  the 
temperature  field  must  be  monitored  on  the  subtime  interval  level.  The  dimen¬ 
sionless  time  step,  Q,  and  the  rate  of  conductivity  change,  (3,  must  both  be 
considered  when  choosing  the  step  size  so  the  solution  does  not  become  unsta¬ 
ble.  For  this  method,  the  term  [ykl?/(m  +  1)]  in  the  exponential  was  consid¬ 
ered  at  its  maximum  possible  value  and  the  time  step  was  adjusted  to  retain 
stability.  This  criteria  was  chosen  so  that 
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A  comparison  of  results  obtained  by  using  a  pure  explicit  method,  a  pure 
implicit  method,  and  the  exponential  method  can  be  found  in  figure  3  and 
table  IV.  Figure  3  shows  the  temperature  field  over  a  slab  cross  section. 

From  this,  it  is  evident  that  the  exponential  and  pure  explicit  methods  give 
very  similar  results.  The  implicit  method  predicted  higher  temperatures  closer 
to  the  slab  surface  and  lower  temperatures  at  the  slab  centerline.  In  table  IV 
the  results  at  the  slab  center  are  shown  for  various  elapsed  times.  As  can  be 
seen,  all  three  methods  agreed  with  each  other  to  within  a  few  percent. 


Unsteady  Heat  Conduction  in  Three  Dimensions 

A  final  application  of  the  exponential  finite  difference  method  to  the 
diffusion  equation  will  be  for  three-dimensional,  unsteady  heat  conduction. 

The  exponential  method,  a  pure  explicit  method,  and  an  implicit  method  (method 
of  Douglas,  ref.  8)  will  be  compared  to  an  exact  solution  for  the  problem  shown 
in  figure  4. 

The  exact  solution  to  the  problem  illustrated  in  figure  4  is  given  in  ref¬ 
erence  9  as 


where  a,  b,  and  c  are  the  widths  of  the  cube  in  the  x-,  y-,  and  z-direc- 
tions,  respectively.  Equation  (33)  will  be  used  to  determine  how  well  the 
numerical  techniques  predict  the  temperature  distribution. 
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The  exponential  finite  difference  technique  will  be  investigated  first. 
The  sequence  to  be  followed  for  determining  the  finite  difference  equation  is 
the  same  as  that  presented  for  the  earlier  cases.  The  step-by-step  procedure 
for  this  three-dimensional  case  consists  of  the  following: 


(1)  Linearize  the  partial  differential  equation 

(2)  Assume  a  product  solution 

(3)  Separate  time  from  spatial  dependence 

(4)  Solve  for  time  dependence 

(5)  Insert  the  appropriate  spatial  finite  differences  into  exponential 
term  that  results  from  step  3 


Based  on  this  procedure  the  three-dimensional  exponential  finite  differ¬ 
ence  equation  can  be  shown  to  be  the  following  (ref.  4): 
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By  using  subtime  intervals,  equation  (34)  becomes 
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where  m  i s  the  number  of  subtime  intervals,  Q  is  the  dimensionless  time 
step,  and  M,  .  .  is  the  dimensionless  drive  number  given  by 
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Equation  (35)  will  be  used  for  all  Interior  node-  in  figure  4.  This  equation, 
as  well  as  those  that  result  from  the  other  analysis,  will  be  modified  along 
the  insulated  boundaries  to  account  for  the  proper  boundary  conditions. 
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The  next  method  to  be  applied  to  this  three-dimensional  case  will  be  the 
pure  explicit  method.  The  finite  difference  equation  for  this  method  is  given 
by  fol lowing  (ref.  8) : 


i  ,j,k 
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1+1.3. k  +  1-1.  j.k  +  1 ,  j+i ,k 


t  n  y  n  -r  n 

i.j-l.k  l.j.k+l  i.j.k- 


where  Q  =  (a  At)/(Ax)  and  Ax  =  Ay  =  Az.  As  shown  in  reference  8,  the 
dimensionless  time  step  Q  must  be 


to  ensure  stability  of  the  method. 


The  last  numerical  technique  that  will  be  applied  is  the  method  of  Douglas 
(ref.  8).  This  method  is  implicit,  and  the  spatial  directions  are  considered 
sequentially  in  the  x-,  y-,  and  z-directions,  respectively.  The  intermediate 
temperatures  U  (found  from  the  x-direction  sweep)  and  V  (found  from  y-direc- 
tion  sweep)  are  used  to  calculate  the  actual  temperature  field  variable  T 
(found  from  z-direction  sweep).  The  equations  that  are  solved  sequentially 
are  presented  as  follows: 
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where  the  finite  difference  operator  in  the  x-direction,  for  example,  would 
be 
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Equations  (39)  to  (41)  must  be  solved  successively  because  the  variable  U  is 
used  in  equation  (40)  to  find  V  and  so  on.  Since  the  method  operates  on  one 


spatial  direction  at  a  time,  the  Thomas  Algorithm  can  be  utilized.  In  the  case 
of  finding  the  U  variable,  the  y  and  z  nodal  values  are  held  constant  for 
the  x-direction  sweep.  This  process  i^s  repeated  until  all  y  and  z  nodal 

values  for  the  x-dlrection  variable  U  _are  calculated.  This  procedure  is 

then  repeated  in  a  similar  way  for  the  V  variable  and  then  finally  for  the 
actual  temperature  field  variable. 

The  results  from  the  three  methods  are  shown  in  table  V.  As  may  be  seen, 
the  exponential  finite  difference  method  gave  more  accurate  predictions  for  the 
nodal  positions  shown.  Also  in  table  V  the  standard  deviation  of  the  diagonal 
values  are  shown.  The  exponential  method  had  a  smaller  standard  deviation  at 
both  elapsed  times  shown  in  table  V. 

In  reference  8  nine  different  methods  to  solve  the  diffusion  equation  in 
three  dimensions  were  investigated.  The  method  of  Douglas  was  the  preferred 
method  because  of  its  accurate  results  and  low  computer  CPU  time.  In  that 
study  the  pure  explicit  method  required  the  lowest  amount  of  CPU  time  with  the 
method  of  Douglas  requiring  approximately  four  times  as  much.  In  the  present 
study  all  three  methods  were  run  on  two  different  mainframe  computers  to  inves¬ 
tigate  how  these  three  methods  compared  in  terms  of  CPU  time.  The  results  are 
shown  in  table  VI.  All  three  methods  were  exercised  for  the  same  number  of 
time  steps.  As  indicated,  the  exponential  method  was  approximately  three  times 
faster  than  the  method  of  Douglas  but  still  slower  than  the  pure  explicit 
method.  From  these  results  it  could  be  concluded  that  the  exponential  method 
would  have  been  chosen  as  the  preferred  method  for  overall  accuracy  and  CPU 
time. 


Viscous  Burger's  Equation 


The  viscous  Burger's  equation  is  given  in  reference  10  as 


9U  „  3U  82U 
31  57  '  "  3? 


(43) 


The  equation  must  be  linearized  first  in  order  to  apply 
Hence,  letting  U  ■  A  =  constant  for  the  nonlinear  term 
equation  give 
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Dividing  by  U  and  evaluating  the  resulting  expression  at  time  n  at  node  i 
result  in 
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The  spatial  and  time  terms  are  now  separated  so  either  side  can  be  set  equal 
to  a  constant  -k 


This  can  be  shown  to  be  equal  to 


Also,  equation  (45)  can  be  shown  to  be  the  following  (ref.  4) 


This  is  used  to  replace  the  exponent  in  equation  (47) 
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Equation  (49)  is  the  exponential  finite  difference  equation  for  the  viscous 
Burger's  equation. 

An  exact  steady-state  solution  to  Burger's  equation  is  available  for  the 
following  conditions: 

U(0,t)  =  U0 
U(L,t)  =  0 

The  steady-state  solution  was  given  as  the  following  (ref.  10): 


U(x)  =  UQU1 


1  -  exp 

V<c  -  O' 

1  +  exp 

U,Re(i  -  l); 

where 

U0L 

Re,  =  —  (5 

L  v 

and  U1  is  the  solution  of  the  following  equation: 

U1  +  1  =  exp  VUlReU 

The  exponential  finite  difference  method  will  be  now  used  to  numerically 
solve  the  previous  problem.  However,  for  the  stated  conditions,  a  problem 
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arises  with  the  portion  of  the  velocity  field  is  Initially  zero.  To  overcome 
this  difficulty,  a  substitution  will  be  used  In  which  a  new  variable  is  defined 
such  that 


U 


Burger's  equation  then  becomes 
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with  the  following  imposed  conditions  if  Ug  =  1 : 


0(0. t)  =  0 

U(L,t)  =  Ug 


(51) 


(52) 


The  same  method  of  separation  of  variables  must  be  performed  on  the  U  vari¬ 
able  in  equation  (51).  The  problem  Is  now  solved  for  the  U  variable  and  the 
substitution  shown  above  Is  then  made  to  find  the  U  variable.  The  exponen¬ 
tial  finite  difference  equation  for  U  can  be  shown  to  be  (ref.  4): 


=  0?  exp 


(53) 


The  results  obtained  by  applying  equation  (53)  and  the  conditions  in  equa¬ 
tions  (52)  are  compared  to  the  steady-state  exact  results  of  equation  (50)  and 
are  shown  in  figure  5.  The  results  from  the  exponential  method  were  nearly 
the  same  as  a  those  from  the  exact  method. 

Another  application  of  Burger's  equation  was  done  to  investigate  the 
effect  of  the  diffusion  term.  The  results  for  the  variation  of  v  over  four 
orders  of  magnitude  are  shown  in  figure  6  for  the  same  instant  in  time.  At 
the  two  lower  v  values,  the  total  range  of  the  field  variable  takes  place 
over  a  small  number  of  nodal  positions.  A  better  approximation  could  be  made 
for  these  cases  by  using  a  finer  grid.  A  comparison  of  the  exponential  and  a 
pure  explicit  finite  differencing  schemes  for  Burger's  equation  are  shown  in 
figure  7.  As  can  be  seen  from  figure  7,  the  number  of  nodes  used  can  have  a 
large  effect  on  the  predicted  velocity  field.  The  pure  explicit  techniques 
can  have  large  oscillations  and  predict  physically  impossible  results.  As  the 
number  of  nodes  are  increased  and  the  time  step  decreased,  the  two  solutions 
give  similar  results. 


Laminar  Boundary  Layer  on  a  Flat  Plate 


The  last  application  to  be  Investigated  will  be  for  the  development  of  a 
laminar  boundary  layer  on  a  flat  plate  (fig.  8).  In  reference  9  the  steady- 
state  formulation  Is  given  In  terms  of  the  following  three  partial  differential 
equations: 


Continuity: 


Momentum: 


Energy: 
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with  the  following  boundary  conditions: 

U(x ,0)  =  0  U(0,y) 
V( x ,0)  =  0  V( x ,  L) 
T<  x ,0)  .  0  T(0,y) 


(54) 


(55) 


(56) 


(57) 


where  v  and  a  are  the  momentum  and  thermal  dlffusivlties,  respectively. 


Equations  (55)  and  (56)  can  be  solved  by  using  the  method  presented  for 
the  viscous  Burger's  equation.  The  only  difference  is  that  the  solution  will 
march  in  the  x-direction  instead  of  time.  The  results  from  the  separation  of 
variables  for  equations  (55)  and  (56)  were  found  to  be  (ref.  4) 
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The  continuity  equation  is  written  as  (ref.  10) 


Equations  (58)  and  (59)  are  first  solved  using  a  spatial  subincrement  as 
was  done  for  the  cases  when  time  was  the  marching  direction  of  the  solution. 
After  this,  the  continuity  equation  <eq.  (60))  Is  solved. 

The  results  of  this  application  are  shown  In  figure  9  for  a  Prandtl  number 
equal  to  0.72.  The  thermal  boundary  layer  was  outside  the  velocity  boundary 
layer,  as  would  be  expected.  The  results  with  the  Prandtl  number  equal  to  0.72 
were  compared  to  the  exact  solution  as  presented  In  reference  9.  A  downstream 
position  was  chosen  and  the  results  are  compared  In  table  VII.  The  exponential 
method  results  were  In  good  agreement  with  the  exact  results. 


CONCLUDING  REMARKS 

In  conclusion,  an  exponential  finite  difference  technique  has  been 
extended  to  other  coordinate  systems  and  expanded  to  model  problems  In  two  and 
three  dimensions.  The  method  has  direct  application  to  linear  partial  differ¬ 
ential  equations  such  as  the  diffusion  equation  and  can  be  extended  to  solve 
nonlinear  equations.  The  method  Is  presented  as  an  alternative  method  for 
solving  a  wide  range  of  engineering  problems. 

The  method  was  applied  to  a  variety  of  heat  conduction  and  fluid  flow 
problems.  It  was  found  that  the  results  predicted  by  the  exponential  finite 
difference  algorithm  for  the  cases  presented  In  this  study  demonstrated  that 

1.  Field  variable  was  predicted  with  a  higher  degree  of  accuracy  than 
other  numerical  techniques  where  exact  solutions  exist. 

2.  The  method  can  be  applied  to  linear  and  nonlinear  partial  differential 
equations  with  dependent  variables  that  can  be  separated. 

3.  When  the  exponential  method  Is  applied  to  the  diffusion  equation,  the 
stability  of  the  method  Is  the  same  as  that  of  pure  explicit  methods,  where 
the  subtime  Interval  step  size  determines  the  stability. 
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TABLE  I.  -  COMPARISON  OF  RESULTS  FOR  DIFFERENT  DIMENSIONLESS  TIME  STEPS  FOR  ONE-DIMENSIONAL 
HEAT  TRANSFER  IN  CYLINDRICAL  COORDINATES  WITH  INFINITE  HEAT  TRANSFER 
COEFFICIENT  AT  THE  SURFACE 

nitial  and  boundary  conditions  are  the  following:  h  -»  so;  T (  r , 0 )  =  1.0:  T ( R ,  t )  =  0.0  for 
t  >  0;  Q  =  (a  At)/(Ar)z;  a  =  1.0  m'vs;  N  =  number  of  nodes;  m  =  number  of  subtime  intervals.] 
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Exponential  finite  difference  results,  °C 
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TABLE  II.  -  FINITE  HEAT  TRANSFER  COEFFICIENT  CYLINDRICAL  COORDINATES 
[T(r.O)  z  1.0,  Tb  =  O,  (2  z  (a  At)/(Ar)2.] 


Time,  Biot 
t ,  modu- 
s  1  us 


Spatial  Exact 

coordi-  analysis 
nate,  (ref.  5), 
r,  °C 


Exponential  finite  difference  results,  °C 
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N  =  11, 
m  =  4, 
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TABLE  III.  -  COMPARISON  OF  EXPONENTIAL  FINITE  DIFFERENCE 
METHOD  IN  ONE-DIMENSIONAL  CYLINDRICAL  COORDINATES 
TO  THE  RESULTS  OF  REFERENCE  6 
[a  =  1.0  ft^/s;  At  =  1.0  s;  (a  At)/Ar^  =  1.0;  N  =  number 
of  nodes;  m  =  number  of  subintervals.] 
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Results  fi 
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TABLE  IV.  -  COMPARISON  OF  EXPONENTIAL,  PURE-EXPLICIT,  AND  IMPLICIT  FINITE 
DIFFERENCE  METHODS  FOR  ONE-DIMENSIONAL,  UNSTEADY  HEAT  TRANSFER  WITH 
TEMPERATURE-VARYING  THERMAL  CONDUCTIVITY 
[Temperature  shown  is  at  center  of  slab;  K(T)  =  1.0  +  B(T);  (3  -  0.01.] 
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33.41929 
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TEMPERATURE 

(B)  LINEAR  RELATIONSHIP  BETWEEN  CONDUCTIVITY  AND  TEMPERATURE. 
FIGURE  2.  -  SKETCHES  SHOWING  PROBLEM  STATEMENT  FOR  TEMPERATURE- 
VARYING  THERMAL  CONDUCTIVITY. 
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FIGURE  3.  -  COMPARISON  OF  METHODS  FOR  TEMPERATURE -VARY  1N6 
CONDUCTIVITY.  SHOWING  TEMPERATURE  FIELD  AT  t  -  0.02  S. 
k(T)  =  kR<1  ♦  pn  WHERE  kR  «  1.0.  P  =  0.01.  T(x.O)  *  100. 
AND  T(O.t)  -  T(L.t>  -  0;  Q  -  1. 


FIGURE  5.  -  COMPARISON  OF  STEADY  STATE  SOLUTIONS  COMPARING 
THE  EXACT  RESULTS  TO  THE  EXPONENTIAL  FINITE  DIFFERENCE 
SOLUTION  U(O.t)  *  1.0;  U(L.I)  »  0.0. 


BOUNDARY  CONDITIONS:  t  >0 
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jT-Cx.y.O.t)  =  0 


Td.y.z.t)  -  T(x.l.z.t)  ■  TIx.y.1,  t)  =  0 

FIGURE  M.  -  BOUNDARY  AND  INITIAL  CONDITIONS  FOR  THREE-DIMENSIONAL 
UNSTEADY  STATE  CONDUCTION  HEAT  TRANSFER. 
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FIGURE  6.  -  EXPONENTAL  FINITE  DIFFERENCE  RESULTS  FOR  VARY¬ 
ING  KINEMATIC  VISCOSITY.  ALL  VELOCITIES  ARE  SHOWN  FOR 
N  =  21.  m  =  8,  At/Ax2  =  M.O.  t  =  1.0  S,  U(O.I)  =  1.0. 
U(L.t)  =  0.0. 
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FIGURE  7.  -  COMPARISON  OF  EXPONENTIAL  AND  PURE  EXPLICIT 
FINITE  DIFFERENCE  METHODS.  ALL  RESULTS  SHOWN  FOR 
V  *  0.01  m2/s.  t  =  1.0  S.  U(O.t)  =  1.0.  AND 
U(L.t)  =  0.0. 
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FIGURE  8.  -BOUNDARY  LAYER  DEVELOPMENT  ALONG  A  COOLED  FLAT  PLATE. 
CONDITIONS:  U(x.O)  =  0.  Vtx.O)  =  0.  VCx.L)  =  0.  UCO.y)  =  1.0. 
T(O.y)  =  1.0. 


DISTANCE  DOWN  THE  FLAT  PLATE.  X,  CM 

FIGURE  9.  -  EXPONENTIAL  FINITE  DIFFERENCE  RESULTS  FOR  BOUNDARY  LAYER  EQUATIONS  WITH  CONDITIONS  U(O.y)  =  1.0, 
U(X.O)  *  0.  V(X.O)  =  0.  VCx.L)  =  0.  T(X.O)  -  0,0,  T(O.y)  -  1.0;  V  =  0.0072  CM2/S;  0  =  0.01  CM2/S. 
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